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We present a fourth order convergent (2 + 1) numerical code to solve the Teukolsky equation 
in the time domain. Our approach is to rewrite the Teukolsky equation as a system of first order 
differential equations. In this way we get a system that has the form of an advection equation. 
This is used in combination with a series expansion of the solution in powers of time. To obtain 
a fourth order scheme we kept terms up to fourth derivative in time and use the advection-like 
system of differential equations to substitute the temporal derivatives by spatial derivatives. A local 
stability study leads to a Courant factor of 1.5 for the nonrotating case. This scheme is used to 
evolve gravitational perturbations in Schwarzschild and Kerr backgrounds. Our numerical method 
proved to be fourth order convergent in r* and 6 directions. The correct power-law tail, ~ l/t"^^^^ , 
for general initial data, and ~ 1/ i^^^"* , for time symmetric data, was found in the simulations where 
the duration in time of the tail phase was long enough. We verified that it is crucial to resolve 
accurately the angular dependence of the mode at late times in order to obtain these values of 
the exponents in the power-law decay. In other cases, when the decay was too fast and round-off 
error was reached before a tail was developed, the quasinormal modes frequencies provided a test 
to determine the validity of our code. 
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I. INTRODUCTION 



In recent years, there has been an increasing interest in 
numerically solving General Relativity's (GR) field equa- 
tions to provide an accurate description of the gravita- 
tional radiation generated in different astrophysical sce- 
narios. This is due to fact that direct measure of grav- 
itational waves will soon be possible with large inter- 
ferometers such as LIGO and LISA. The astrophysical 
events that produce gravitational waves are among the 
most energetic phenomena ever seen. The best candi- 
date for such observations are the collisions of binary 
black hole systems. In order to detect this gravitational 
radiation, accurate templates are needed, and this im- 
plies to solve the non-linear GR equations. This proved 
to be a very challenging task. Although in the last or- 
bital stages (merger) of a binary black hole collision it 
is necessary to solve Einstein's equations (full numerical 
approach), the very last part of the coalescence (close 
limit) the system can be considered as a single distorted 
black hole and can be treated with perturbation theory. 
Comparisons between full numerical simulations and per- 
turbative methods show a surprising agreement. This has 
encouraged people to go beyond matching close limit and 
full numerical simulations in tandem to produce simula- 
tions that neither of each technique alone was able to 
do. The general method of coupling full numerical and 
approximate techniques is the main development of the 
'Lazarus project' [1|1I1I1I3- 
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Moreover, there are physical scenarios in which the 
metric departure from the known static black hole solu- 
tions is always small outside the horizon. Some examples 
are a particle (compact object) orbiting a black hole, the 
propagation of gravitational waves and accretion disks 
near black holes. The equations evolved in the pertur- 
bative method are linear which is an advantage over the 
non-linear set of ten coupled general relativistic equa- 
tions. 

One application of perturbation theory is the propa- 
gation of waves in a curved spacetime. In this work, we 
will focus our study in the late time behavior of a grav- 
itational wave in such curved spacetimes. Although the 
problem is not new, according to Andersson l^j there are 
still some aspects for which there is no definite answer or 
no answers at all; like the role played by highly damped 
modes and the intermediate behavior of the power-law 
tails. Power-law tails is the name given to the last phase 
in the propagation of a wave in a curved background, 
as we will see below. The general features of the evolu- 
tion of such a test field in the proximity of a black hole, 
as seen by a distant observer, can be divided in three 
stages: i) Radiation emitted directly by the perturbation 
source. It depends on the form of the initial field (initial 
data) . ii) Quasinormal ringing. It depends on the param- 
eters of the black hole. They are exponentially damped 
oscillations and carry part of the gravitational radiated 
energy in astrophysical processes like gravitational col- 
lapses. Quasinormal modes are characterized by complex 
frequencies a. Such modes can be represented as e**^*, the 
real part of cr is the oscillation frequency and the imagi- 
nary part is the exponential damping, iii) Power-law tail. 
The field decays with time according the a power-law at 
very late times. 
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Most of the work on the late time radiative fall-off 
has been done in spherically symmetric spacetime. In 
this case the solution admits a decomposition in spherical 
harmonics. The field equations are reduced to a wave 
equation with a effective potential. For the Schwarzschild 
case this wave equation is 
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^i{t,r)^0 (1.1) 



where r* is the tortoise coordinate 



and 



r* =r + 2M\n(r/2M -I), (1.2) 
l{l + 1) 2M' 



(1,3) 



being M the mass of the black hole. 

In 1972, Price 't'I treated Eq. (|1.1|) as a perturbative 
expansion in powers of M . He showed that the perturba- 
tion decays in time according to ^-(2^+3) foj- 1 '3> r ^ M. 
This situation occurs at a finite value of r while t — > 00, 
which is called timelike infinity. 

Quasinormal modes and power-law tails have been 
studied by Leaver who analyzed the problem in the 
frequency domain. He found the correct late time be- 
havior with a low-frequency approximation. Quasinor- 
mal modes can be considered as the "pure tones" of the 
black hole. Once excited, their damping and oscillation 
frequency depend only on the parameters of the black 
hole. Ching et. al. argue that the late time decay of 
the field can be seen as a scattering due to the spacetime 
curvature. This implies that the power-law behavior de- 
pends only on the asymptotic conditions of the metric. 
In a recent work, Poisson 10] found that in a weakly 
curved spacetime the late time dynamics is insensitive 
to the non spherical aspects of the metric, being entirely 
determined by the spacetime total mass. 

In Ref. [ij it is found that the power-law tails are 
not an universal phenomenon as thought, but the fall-off 
power depend on the initial profile of the fields. While 
for the generic data with <i>£|*=" ^ and dt^e\^=° 7^ the 
predicted decay goes like l/t^^"*"^, when you start with 
an initially static field, i.e. ^ and at$£|*=° = 

the predicted decay goes like l/i^^^"*. 

Most of the work with the Teukolsky equation has 
been performed in the frequency domain (quasinormal 
modes, wave scattering, motion of test particles). Per- 
turbation on the frequency domain can be reduced ana- 
lytically to solve ordinary differential equations and can 
lead to a better understanding of the physics involved in 
the phenomena. Such information is much more difficult 
to obtain from purely numerical computations. However, 
in more complete treatments, the number of frequencies 
that one needs to consider is orders of magnitude larger 
than the number of points needed to resolve the di- 
rection. Furthermore, the study of quasinormal modes 
would require higher resolution near w = to resolve the 



tails. The resolution of the quasinormal modes is also 
sensitive to the spacing in frequencies. These are the ar- 
gument presented by Krivan et. al. [3| in favor of a 
numerical treatment of perturbations using the Teukol- 
sky equation in the time domain. Another motivation to 
work in the time domain is that when one tries to find a 
solution of the radial Teukolsky equation in the frequency 
domain with a source term that extends to infinity, the 
result is divergent. This means that the Teukolsky equa- 
tion needs to be regularized when sources are present 

[Hill. 

The difficulty with numerical integrations of the 
Teukolsky equation is the linear term in s on the first 
time derivative. Depending on the relative sign between 
this term and the second time derivative, it can act as 
a damping or anti-damping term. As described in jl2|. 
good time evolutions are achieved by writing the Teukol- 
sky equation as a set first order differential equations 
system. This and other issues concerning the implemen- 
tation of a fourth order algorithm are discussed in the 
next chapter. Motivations for a fourth order convergent 
algorithm are mainly, that it can reproduce the same 
results of a second order convergent code with less res- 
olution. This makes the former to run faster than the 
later. Although in a fourth order scheme, more intense 
computation is needed, there is some gain in speed when 
equivalent resolutions, i.e. resolution that produce the 
same error in the solution, are used. On the other hand, 
if equal resolutions are used, a fourth order method will 
yield more accurate solutions. The price we have to pay 
in this case is that the fourth order method will take 
longer. This issues enter in consideration in gravitational 
wave detection, because a large number of templates need 
to be generated. The gain in accuracy can also be used 
in second order perturbation theory, since higher order 
derivatives of the field are needed to build up the eflFective 
source term |l5j| . This may have important applications 
in the 'Lazarus approach' 0| and the radiation reaction 
problem of a particle orbiting a black hole. 

Using the Kinnersley null tetrad: 
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[(r2+a2)/A,l,0,a/A] 
= [r2+a2,-A,0,a]/(2E) 
^ = [ia sinO, 0, 1, i sin 0]/[V2{r + iai 



(1.4) 



where T. ^ + a"^ cos^ 9 and A = - 2Mr + The 
Newman-Penrose equations written in Boyer-Lindquist 
coordinates lead to the Teukolsky equation 
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The gravitational perturbations are recovered setting s = 
±2, s is a parameter called spin weight. Here the quantity 
of interest is ^' = p'^'ipf with p = — l/(r — ia cos 9) and 
s = -2. For is s = 2, = -0^. 

In vacuum, it is well know that 1)1. 5|1 can be sep- 
arated in the frequency domain, by taking 5* = 
e-i^t^',"^4>s{9)R{r). When s = 0, the functions S{e) 
are the spheroidal functions. When aoj = these 
eigenf unctions are the spin weighted spherical harmon- 
ics sYi'^ie^cj)) =s S'["(6l)e™*. The general solution can 
be represented as 
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In most astrophysical applications, we are interested in 
computing solutions that represent gravitational radia- 
tion at infinity. This information is carried by . This, 
together with the good asymptotic behavior of the solu- 
tions are the motivations of choosing the value s = — 2 to 
perform the time evolutions. Knowing the value of 
allows us to calculate the outgoing energy flux per unit 
time asM^s! 
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where dil = sin 9 dO d(j) and u = t — r. 

In the next section we explicitate the numerical tech- 
niques used to solve the Teukolsky equation starting from 
a review of the second order formalism developed in 
Ref. We then derive the forth order accurate in 

time and space algorithm, as a natural generalization of 
the Lax-Wendorff algorithm. It is followed by a stability 
analysis that find a Courant factor of 1.5 between the 
time and spatial coordinate r* that ensures stable evo- 
lution. We finish Section II with a description of the 
radiative boundary conditions imposed to the field in the 
two radial boundaries and the angular boundary condi- 
tions imposed by the symmetry of the solution. Section 
III deals with the applications of the fourth order code 
developed. As initial data we consider an outgoing Gaus- 
sian pulse located in the far region. The fourth order 
convergence in both spatial variables in verified in Sec- 
tion III.B. We this we compute the quasinormal modes 
and power-law tails and an ultimately test of the code. 
Finally we discuss some of the consequences of the com- 
puted decay powers on the light of some discussion in the 
literature. 



the time domain by rewriting it as a set a first order 
partial differential equations. The aim of their work is 
the evolution of gravitational perturbations of the Kerr 
spacetime. This is done for the case s = — 2. This choice 
is due to the fact that for this particular value of spin, 
the solutions are bounded for both, ingoing and outgoing 
waves near the horizon (r* — > — oo) and at infinity (?'* — > 
oo). The asymptotic behavior of the solutions of (|1.5|l . for 
a given spin weight s, representing ingoing and outgoing 
waves is [l3, [lal 



lim 

r* — >-t-oo 



lim jvps 

r* — > — oo 



2^j.2s+i outgoing 
1/r for ingoing 



1 for outgoing 
A^'' for ingoing 



(2.1) 



(2.2) 



The procedure of solving the Teukolsky equation is ini- 
tiated by introducing the ansatz 

*(t,r*,6l,0) = ^r3e™^$(t,r*,6') (2.3) 



where r* is the Kerr tortoise coordinate, defined as 
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and (b is the Kerr azimuthal coordinate, defined as 



dcf) — d(j) + — dr. 
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Now they introduce an auxiliary field 11 defined as 
U = dt^ + bdr'<i>, (2.. 

where 



b = 



r^+a^ 



(2.9) 



II. NUMERICAL TECHNIQUES 
A. Review 

In this work we will follow the approach used by Krivan 
et. al. They integrate the Teukolsky equation in 



E^ = (r^ + a^)^-a^Asin 



(2.10) 



This decomposes the second order differential Teukolsky 
equation (|1.5|) into a system of first order differential 
equations. The resulting equation has the form 



dtu + Mdr'U + Lu + Au = 



(2.11) 
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where u = (<I>ij, $/, 11^;;, 11/)^ is a column vector whose 
components are the real and imaginary (i?, I) parts of 
the fields $ and 11. The coefficients of the derivatives of 
(|1.5|) are rearranged as the elements of matrices M and 
A, given by 



(2.12) 
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131 132 
—132 031 



033 034 
—034 033 



(2.13) 



The remaining matrix L contains the derivatives of the 
fields with respect to the polar coordinate 9, whose com- 
ponents are 







/31 

^31 



(2.14) 



The deduction of substitution H2.8|l and the value of each 
element of the above matrices are given later in this chap- 
ter. 

They proceed to solve the first-order system (I2.11|) us- 
ing the Lax-Wendroff method jT^. For this, l|2.11|l is 
written as 



dtu + Ddr'U = S, 

-b, —b) and S = 



(2.15) 



where D = diag(6, 6, 
Lu. 

To solve the above equation in Ref. |0 it was used a 
grid with 8000 points for r* and 32 points for 6. The com- 
putational domain was -lOOM < r* < 500M and < 
9j < TT, with a Courant condition of St < min( 6r* , 5 SO). 

Boundary conditions have been imposed as follows: 
$ = n = at the horizon and outer boundary. Along the 
axis, <I> = or 9^$ = for m, (the azimuthal number) 
odd or even, respectively. 

With this settings Ref. ^3 code showed stability of 
the order lOOOM of evolution time. It was second order 
convergent for times < 50Af , with a convergence rate 
higher than 1.3 for later times. 

In the following sections we will concentrate in general- 
izing the numerical algorithm to accomplish a fourth or- 
der convergent numerical evolution, using the first-order 
formulation of the Teukolsky equation H2.11II as starting 
point. 



B. Rewriting Teukolsky equation 

1. Separating dependence 

The Teukolsky equation is separable in the azimuthal 
variable for the general case of s 7^ and a 7^ 0. Fur- 



thermore, we will change the normal Boyer-Lindquist co- 
ordinates r and by r* and 0, respectively. This has 
the advantage that r approaches asymptotically to r+ as 
r* goes to minus infinity, so the inner boundary of the 
computational domain is approximated very close to the 
horizon, but still outside the black hole. 

The variable </> is used to improve the behavior of the 
coordinates near the horizon. This is a manifestation of 
the frame dragging effect of a rotating black hole • 

Now we make use use of the ansatz (|2.3f) used by Kri- 
van et. al., where they also include a factor r^. This is 
done to eliminate the increasing behavior of the solutions 
at infinity, according to (|2.1|) . So we substitute H2.3|l into 
(|1.5|) , setting the source term T = to get vacuum space 

solutions. After dropping a global r^e™"^ factor we get 



2^2 



(r" + a ^) 



- sin^ ( 



dtt^ 



Ms{r^ - a^) - ^-sA 



ia(sA cos 6 + 2Mmr) 



1 



ft* 



A 

A 



+ — r KSr^ + 6a^)A - 2rs{r^ + a^){M - r) 



)] dr-^ + dee^ + cotede'^ 



6Af(s + l)-r(7s + 6) 



BA-" - rA 
+ r(s cot 6 + m esc 9)^ 
- 2iamr [2rs{M - r) - 3A] | = 

We can bring this equation into the form 



(2.16) 



+Ceedgg^ + Cede^ + Cso* = (2.17) 

where the C's represent the coefficients of the derivatives. 
We simple multiply by — A/S^. The results are: 



Ct = 2s 



2ia 



2mMr + s A cos ( 

S2 



M{a^ - r^) + rA 
rs{M - r)(r2 + a^) - {ia"^ + 4r2)A 



rS2 



—2iam- 



S2 



Co 



■ cot( 



(2.18) 
(2.19) 

(2.20) 
(2.21) 
(2.22) 



Cso — 2iam 



2rs{M - 



-3A 



rE2 



6Mr{s + 1) - r2(7s + 6) - 6A r2(s cot 6* -h mcsc 61)2 



r2E2 



(2.23) 
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With the Teukolsky equation written as (|2.17() , we may 
now proceed to transform it into an equivalent first-order 
differential equations system. 

2. From second to first-order differential equations 

A common technique to numerically solve a second- 
order differential equation is to rewrite it as a set of 
first-order differential equations and apply the appropri- 
ate methods. For cases such as the simple wave equation, 
getting such a first-order system is straightforward. The 
equation reads 

dttu = vld^.u, (2.24) 

with Vp the velocity of propagation and u = u(t,x). It 
can be written as 

duu - / (2.25) 
/ = vld^^u. (2.26) 

Here / is a function chosen in such a way that second 
derivatives are eliminated. In this case f = Vp dtxV will 
do the job, with v = v(t, x), and we have 

dttu = VpdtxV (2.27) 
VpdtxV = VpdxxU- (2.28) 

Now we factor out the time derivative in the first equation 
and the spatial derivative in the second one and we get 
a system of first-order equations 

dtu — VpdxV (2.29) 
dtv — VpdxU (2.30) 

or in matrix form 



u 




■ 


Vp 




u 


V 


t 


Vp 







V 



(2.31) 



where the subscripts denote derivatives. Notice that it 
has the form of the advection equation for a system of 
equations 

dtU = Yp dxu (2.32) 

where u = {u,v)'^ is a column vector. 

Trying to carry out the above procedure for the Teukol- 
sky equation is more complicated and probably will not 
yield results, because it contains cross derivatives. Our 
approach here will be different. Recall Eq. H2.17|l . in 
which we set the coefficient of the second time deriva- 
tive equal to one, dividing the equation by precisely this 
coefficient (the subscripts in the C's do not mean differ- 
entiation with respect to the variables). 

We want to write H2.17|l as an advection equation with 
a source term: 



" $ ■ 








■ $ ■ 




711 712 






n 


-1- 

t 


P21 P22 




n 


-1- 

r* 


721 722 




n 



= 0. 



(2.33) 



where none of the /3's and 7's depend on t. Expanding 
the matrix products we get two equations 

9t$-K/3iia^.$-K/3i29r*n-K7ii$ + 7i2n = o, 

(2.34) 

dtIl + /32ldr''S>+l322dr'Il + j2l'S>+l22ll = 0. 

(2.35) 

Now let Pi2 = 7ii = and 712 ~ —1. By doing this 
we obtain an expression for 11 that depends only on the 
derivatives of $ with respect to t and r* (and some func- 
tion /3ii): 

n = 9t$-F/3iia,.$, (2.36) 

which is similar to (|2.8|l . but Pu has not yet been speci- 
fied. We will see that with this definition, 11 can be easily 
substituted and eliminated in the second equation H2.35|l . 
The derivatives of 11 are 

dtU = du^ + f3iidtr'^ (2.37) 
dr.U = a^.t$-H(5^./3ii)(9,..$)-l-/3ii5r.r-$(2.38) 

Substituting ifT^ . ifT^ and into ifT^ and re- 

arranging terms we find 

dtt^ - I22dt^ - iPn + /322)9r.t$ + l322Plldr' 

-(/321 - I322dr'f3ii - 722/3ii)5..$ - 721$ = (2.39) 

To find the value of the /3's and 7's we equate the coeffi- 
cients of the derivatives of $ in (|2.39l) with those in H2.17|l . 
From the coefficient of dr't^ we find that /322 = —Pii- 
Combining this with the coefficient of dr*r'^ we see that 

Pll — \/—Cr-r' ■ 

- P22 = Pii = ^^^^ = b (2.40) 

This is exactly the definition of b in (|2.9|l . It is easy 
to show that the remaining equations yield the following 
results: 

722 = Ct (2.41) 
721 = Cso + ^ai (2.42) 
^21 = C'r' +bdr*b-Ctb. (2.43) 

Here ^31 is defined according to H2.14|l as the 6'-derivative 
operator: 

hi^Cggdee + Cedg, (2.44) 

in the sense that it can be "factored" as ^31$ and added 
to Cso'i'- Of course, an expression as H2.42|l is not mathe- 
matically rigorous. It is rather a way to express that the 
equation's angular dependence is going to be added to 
the source term. Furthermore, in a numerical implemen- 
tation we don't compute the value of 721, but the value 
of 721 
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3. Splitting real and imaginary parts 

Because Teukolsky equation involves complex coeffi- 
cients it is necessary to treat the real and imaginary parts 
of the solution. Let's define four functions 11^ 
and 11/, such that 



Summarizing, the coefficients used in this work, using 
the ansatz (|2.3|l are 



P^i = 2 



rs{M - r){r^ + a^) - (Sa^ + Ar'^)A 
2bs ^ +bdr>b (2.53) 



$ = $ 

n = Hr 



(2.45) 
(2.46) 



S2 



and substitute them into H2.34|l and H2.35|l . After col- 
lecting real and imaginary parts and equating both to 
zero, we obtain a set of four equations. As a shorthand 
to denote derivatives, we use a dot for dt and a prime for 



P21 = —'2am 
A 



"2 ' "2 2mMr + sAcos( 



4- 



26a- 



S2 



2 , 6Mr{s + 1) - r^{7s + 6) - 6A 

r'^{s cot 9 + TOCSC^)^^ 

A A 

:dee — cot 9de 



(2.54) 







$7 + 6$j - H/ = 



(2.47) 



(2.48) 



2am 



2rs(M - r) - 3A 



,2 



Cf = 2, 



nK+/32''i*fl-/32i*/-^n',,+7fi$fl-72\$/+c/^nfl-c/n/ = o 

(2.49) 
(2.50) 

We can finally arrange these equations in matrix form as 
follows: 



= 2a 



2mMr + sA cos 9 
S2 



(2.55) 
(2.56) 
(2.57) 
(2.58) 



$7? 



+ 



6 

6 

/3fi -b 

/32'i -6 









-1 







' $7? ' 











-1 




$7 


721 


-7I1 




-cl 






7I1 


721 


ci 


_ 




L n7 J 









$7 








. n7 . 




0. 



C. Derivation of the 4th order algorithm 

Let's assume that u{t, x) is a continuously differen- 
tiable function and that its derivatives in both t and x 
exist up to order four in some given interval. This func- 
tion satisfies a differential equation of the form 



Ut = VUx 



(2.59) 



(2.51) 



Comparing this last equation with 12.11|l and making 
use of our definitions of (32i and 721, equations (|2.43() and 
(|2.42|) . respectively; we see that our derivation yields the 
same structure for Teukolsky equation as that derived 
by Krivan et. al. It is worth to comment here that 
although the structure of the equation is the same, the 
coefficients shown above do not agree completely with 
those of Ref. . The coefficients that do not agree are 
Cso and the real part of Cr* , which correspond to what 
Ref. ^2 calls C2,C5 and cg. The coefhcients reported in 
[T^ correspond exactly to the case in which the ansatz 
for the solution is 



where u is a constant and the subscripts represent first 
derivatives with respect to t and x respectively. The 
Taylor expansion in t for u{t, x) is 

St"^ 6t^ Sf^ 
u{t+ St,x) = u + utSt + utt-^ + uttt-^+Utttt^-\ ■ 

■ (2.60) 

Here u and its derivatives are evaluated at some given 
point {tojXo). Now we can use eq. (j2.59|l to replace the 
time derivatives by spatial derivatives, using the fact that 
Utt = v'^Uxx, uttt — v^Uxxx and so on. Thus we have 



u{t+6t,x) = u + VUx 6t + v'^Ux 



2! 



^{t,r\9, 



3'""^<^{t,r*,9). 



(2.52) 



Notice that this equation does not contain the function 
used in (|2.3|l . Another correction that needs to be done 
in order to recover Teukolsky equation is to set 031 = C5 
and 034 = — C3, according to the definitions presented in 



3 St^ 4 / ^ 

V Uxxx~^^^ \~ V Uxxxx~^^^ ^ ■ ■ {2.61) 

If we truncate this series, taking terms up to St^ we will 
obtain the Lax-Wendroff scheme, which is second order 
accurate in space and time. Introducing the usual dis- 
crete notation [/" for u{t„, Xj) and using a second order 
accurate approximation for x derivatives we have [l9l | 



n+l 



^3 



-a 



u: 



(2.62) 



where a = v 5t/ Sx. 

In order to have a fourth order accurate scheme, we 
truncate the Taylor expansion of x) including terms 
up to 5t* and using the differential equation H2.59|l to 
replace time derivatives by space derivatives. Further- 
more, we use a fourth order finite difference scheme (see 
appendix ^) to approximate spatial derivatives. Pro- 
ceeding in this way we have 



2 



12 



a 



24 

144 ^^^-3 

-i2c/;v2 + t/jV3) 



13 13 iJjVi 



12 [/"_2 + 39 C/"_i - 56 [/" -f 39 [/"+! 




k5x 



FIG. 1: 151^ vs. kSx for several values of a. 



(2.63) 



Notice how the dependence on the discretization param- 
eters St and Sx appears as powers of a. 



|^(fc)p = l+gtt^ (cos fc - 5) sin' 



--a"* (4 cos kSx - 17) sin* 



k Sx 
k Sx 



D. Stability analysis 

We apply now a von Neumann stability analysis to the 
previous scheme. This analysis is local, which means that 
we assume that the coefficients of the finite difference 
equation vary slowly in space and time such that they 
can be considered to be constant. In our case, these 
coefficients do not depend on time. We say that the 
method is stable if the scheme is stable for every constant 
value of the coefficients in their range |23]. The idea is 
to expand the solution of the difference equation in its 
eigenmodes e*'"'-' ''^, where fc is a real wave number. The 
time dependence of these modes is a succession of powers 
of some complex number ^(fc) , called amplification factor. 
With this, we say that the difference equation is stable if 
< I, for a given value of fc. The eigenmodes of the 
difference equation H2.63|) can be written as 0| 



^3 



n Akj 5x 



(2.64) 



where ^ = ^(fc) is a complex quantity and fc is a real 
wave number. Substituting H2.64|l into H2.63|l we get a 
first grade polynomial in ^. After some algebra we get 

^(fc) = 1 + ia^(cosfc(5a; - 7)sin^ - 



-a (cos fc Ox — 4) sm — — 



a(8 sin kSx — sin 2fc Sx)~ 



— Q;'^(8sin2fci5a; — 13sinfci5x — sm3kSx) 



(2.65) 



— a® (133 cos fc — 34 cos 2 fc Sx+ 
54 



3 cos 3 fc (5a; — 150) sin* 



fc Sx 



4 o , ^2 ■ p. k Sx 

— a (cos fc ox — 4) sm . 

81 ^ ^2 



(2.66) 



This is a periodic function with a period of 27r. For 
the scheme to be stable it has to satisfy the stability 
condition 



1^1 <1- 



(2.67) 



Although we can find analytic solutions for this equation, 
just to present them here would occupy several pages. 
Rather than analytic solutions, we are interested in some 
interval of values for a, such that our fourth order method 
is stable. The behavior of as a function of fc is shown 
in Fig. n for some values of a. With certainty, we can 
conclude that ^ is less than 1 for a < 1.5. This implies 
that the Courant condition for this case is 



V St < 1.5 Sx 



(2.68) 



This value is 50% greater than that required by the 
Lax-Wendroff method and many others. This can be un- 
derstood by expanding the amplification factor ^(fc) in 
series of fc Sx. In practice Sx must be small enough to 
correctly approximate the continuous differential equa- 
tion. This means that for modes corresponding to small 
values of fc we can expand (|2.66f) in a power series of fc Sx. 
That is 



+ 0((fc(5x)^). (2.69) 



\i\^ = i-[- + -]ikSxr 



8 



It is interesting to compare this with the corresponding 
series expansion of for the Lax and Lax-Wendroff 
methods. For the Lax method we find [T^ 



|eP = l-(l-a')(fc<5x)2 + 



and for the Lax-Wendroff method we have 



(2.70) 



(2.71) 



We can see that our fourth order method has a sixth or- 
der dependence of k 6x, which means that mode damping 
effects become relevant for much higher values of k, mak- 
ing this method more accurate. The generalization of the 
scheme to two dimensions is illustrated in the last sub- 
section. In this case the situation is more complicated 
because we are taking first and second derivatives in the 
6 direction. We were not able to find an analytic Courant 
factor that took into account the 9 direction. From nu- 
merical experiments we verify that the Courant condition 
used in [l^] was reliable in our case too. Thus, as a rule 
of thumbs, we always kept St = min( 6r* , 5 69). 



E. Boundary conditions 

1. Radial boundary conditions 

We use Sommerfeld boundary conditions at radial in- 
finity and at the event horizon. When one uses the tor- 
toise coordinate r* the event horizon is reached when 
r* — > — oo. In practice, it turns that setting r* = 
— 50M is a good approximation. For this value we have 
\r — 2M\ « 10~^^. At the inner boundary, the condition 
is that of an ingoing wave 

-^t,r*,9)^—'i>{t,r*,9). (2.72) 

At the outer boundary, the appropriate condition is that 
of an outgoing wave 



(2.73) 



In order to make this conditions compatible with our 
fourth order integration scheme, we take higher deriva- 
tives of H2.72|l and H2.73|l . The idea is again, to substi- 
tute the time derivatives of the Taylor expansion (|2.t)U|) 
by means of the boundary conditions above. The results 
are summarized in table 

The implementation of boundary conditions is 
straightforward when a second order scheme is employed. 
This is due to the fact that we only need up to second 
order spatial derivatives and its stencil demands only 3 
points[30|. On the other hand, the case of a fourth or- 
der accurate expression for spatial derivatives needs two 
more points for the first derivative and seven points for 
the fourth derivative, see appendix ^ So, in this case. 



llllitil UUUllLlcLl y / — ' min 


vyULtli UULlllLlcLl _y / — 'max 








al* = a^.-i- 











TABLE I: Boundary conditions for the radial direction. 



Off-centered derivatives 


used at points 




1, iVr - 1 




1, 2, Nr-1, Nr-2 



TABLE II: Points for which off-centered derivatives in the r* 
direction are used (point labeling: A^,- + 1 points from to 

Nr). 



the way we implemented the radial boundary conditions 
is to use off-centered expressions to compute the spatial 
derivatives, when needed. Assuming a computational 
grid of A^r + 1 points in the r* direction, labeling the 
points from to A^,-, table ITU shows the points for which 
off-centered spatial derivatives are used. 

As to the auxiliary field 11, defined in 1)2.8(1 . its bound- 
ary condition follows directly from its definition and from 
((2.72(1 and ((2.73() . Once we know the boundary values for 

the value for 11 at the inner boundary is 



n 



(6-1-1)9,-$ 



(2.74) 
(2.75) 



where we have used ((2.7211 and b is given by ((2.9(1 . In a 
similar way, the outer boundary condition is 



(2.76) 



n= (6-1)9,.$ 



2. Angular boundary conditions 

These are imposed along the rotation axis, i.e. at 9 = 
and 9 = TT. The boundary condition depends on the 
particular azimuthal mode m (see equation ((2.3(1 ) chosen 
for the evolution, it can be stated as 

$ = for m = ±1, ±3, ±5... (2.77) 
5e$ = for m = 0,±2,±4... (2.78) 

This conditions come directly from the behavior of the 
solution in the 9 direction. At the same time it is pre- 
cisely this behavior what we use in order to implement 
the appropriate boundary conditions. The solutions for 
which m is even, have even parity about both, 9 = and 
9 = n. On the other hand, the solutions with odd m, 
have odd parity about 9 = and 9 = tt, i.e. 



$(i,r* 
$(^,^*,7^^ 



$(<,r*, 
$(t,r*,7r - 



for m = 0,±2,... 

(2.79) 
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FIG. 2: White circles: ghost zones, black circles: normal grid 
points 



for TO = ±1, ±3, ... 

(2.80) 



To take advantage of this property we need to use a stag- 
gered grid in the 9 direction. By doing this we also avoid 
the inherent difficulties of evaluating expressions in which 
cot 9 is present, (like the last term of Teukolsky equation 
()1.5|l ^ since this function is not finite neither at = 
nor 6 = TT. In a staggered grid, the values for 6* = 
and 9 = n are always located exactly between two grid 
points. The points to the left (right) of 9 = {9 — n) 
are considered as "ghost zones" , because are used just to 
implement the boundary conditions. In our fourth order 
method, we need four ghost points. Two before the first 
point immediately after 9 — and two more after the 
point immediately before 6* = tt, as shown in Fig. 12 Our 
grid has Ng + 3 points in the 9 direction, the first two 
and last two points are ghost zones. In this way we can 
always use a centered formula to compute the derivatives 
in the 9 direction. Notice that if we are using a second 
order accurate approximation, we only need one ghost 
point at each end of the grid. The values of the ghost 
zones are updated according to H2.79|l and (|2.8U|) . 

These properties of the solutions are a direct conse- 
quence of the spherical harmonics behavior (when the 
spin parameter s = 0) and the spin weighted spherical 
harmonics (when s — —2). 



Here we have used the fact that none of the coefficient 
of the Teukolsky equation are time dependent and that 
partial derivatives commute. 

The four time derivatives above could be computed us- 
ing just the finite differences formula for the first deriva- 
tive in r* . Once we know 9tU we can numerically substi- 
tute this result to get the second time derivative &lw and 
so on. The problem with this procedure is that because 
of the exclusive use of the first r* derivative formula; at 
the end, we will not have a fourth r* derivative with the 
accuracy shown in appendix ^ It still will be fourth or- 
der accurate but using the first derivative formula four 
times will propagate more error than that of the fourth 
order accurate finite differences formula. The same holds 
for the other derivatives. The approach we took was to 
compute all time derivatives directly from the coefficients 
of the evolution equation, and the value of the fields at 
every time step. Of course, this implies much more larger 
expressions to compute time derivatives, because now we 
have to algebraically substitute one time derivative into 
the other. Carrying out such a substitutions we find: 



dtu = 
dfu = 



dfvL 



Bu' + Su (2.86) 
B [B'u' + (Su)' + Bu"]-f S9tu (2.87) 

B I (S 9tu)' + B'^ u' B' [(S u)' 3 B u"] 
+ B u' B" + (S u)" + B u(3) I + S a^^u (2.88) 
B {(S52u)' + B"'u' + B'' ((Su)' + 7Bu") + 
B' 



S dtvi)' + B (^4 u' B" + 3 (S u)" + 6 B u^^^ 
B [(S u)' B" + (S dtu)" + B (4 B" u" 
u'B(3) + (Su)(3)+Bu(4))]} 

(2.89) 



F. Notes on implementation 

We mention here some implementation details of equa- 
tion H2.11|l . which for convenience will be written as|3lj 

5tu = -mdr' u - (L + A)u. (2.81) 

As stated above, the main idea has been always to 
substitute the time derivatives in the Taylor expansion 
(|2.t)U|l by spatial derivatives using our differential equa- 
tion ifTHT)) . Calling B = -M and S = -(L A), all 
time derivatives needed are: 

dtVL = Ba^.u + Su (2.82) 

d^vL = mr-{dtVL) + S{dtVL) (2.83) 

d'^tVL = Ba^.(5» -l-S(a» (2.84) 

d'^vL = mr'{d^u) + S{d'^^u). (2.85) 



To clarify the equations, we have used "primes" to denote 
differentiation with respect to r* . 

There is one final issue, worth mentioning here. No- 
tice that there are still some products in which the time 
derivatives of u and u itself appear explicitly. All of 
these products involve multiplication with S. Neither 
the time derivatives nor u are algebraically substituted 
in these products because S = — (L -I- A) and L contains 
the 9 derivatives operator. Instead of expanding further 
derivatives, we chose to numerically substitute the time 
derivatives of u and u itself into these products. We use 
the term "numerically substitute" in the sense that each 
time derivative is calculated and stored in the memory 
of the computer, further computation makes use of the 
stored values. 

Such a procedure gives good results as shown in the 
next section. 



i-modc 

= 1 
= 2 
= 3 
= 4 



constant 
cos 6 
3 cos^ e - 1 



5 cos 



3 cos 6 



35 cos" 61-30 cos^ 9 + 3 



2-'; 



(5 + 7 cos 29) sin^e 



TABLE III: 9 dependence of Y° and -2Y° without the nor- 
mahzation constant. 



III. RESULTS 

A. Initial data 

In all runs, initial data with compact support was used. 
The function used was a Gaussian bell centered at r* = 
75M, in the r* direction and some I, m mode dependence 
in the 9 direction. Thus, for t = we have 
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FIG. 3: Convergence test in r* at 6* = 7r/2 and t = 200M. 



$(0, 



= g-('"*-75)^/i00Q^^j-^-j (3 1) This means that when plotted together, p^.(\E'/i 



where Qim{0) represents the 9 dependence as spherical 
, (j)) or the spin weighted spherical har- 
); for s — Q OY s — respectively. Ta- 



harmonics Y, 



.Y, 



^medium) and medium - "^coarse HlUSt lie OU top of each 

other and that will indicate that our numerical scheme is 
fourth order convergent. The same applies in the case of 



ble mil shows the 9 dependence of the first four spherical 
harmonics and spin weighted spherical harmonics. 



B. Fourth order convergence 

Convergence was tested in both r* and 9 directions. 
The method to assess convergence was to compare three 
runs for the same initial data but different resolutions. 
If we want to measure convergence in the r* direction, 
we keep 9 resolution fixed; while we vary r* resolution. 
The same holds in the case of assessing 9 convergence. 
The way of varying resolutions is such that they keep the 
same ratio. If we call this resolutions fine, medium and 
coarse] they satisfy: 



Starting with the r* direction, Figs. 13151 shows fourth 
order convergence. The simulation parameters are (in 
all runs the black hole mass is taken as M = 1 and the 
Courant factor is 0.5): 



domain: 

grid size [Nr* x Ne). 
resolution: 
physical parameters: 
initial data: 



-50M < r* < 950M 
1000 X 8 

5r* ^ {1,0.5,0.25}M, 59 = vr/S 
a = 0,Z = 2,m = 0, s = -2 
ingoing Gaussian pulse 



5r: 



Sr* 



mediu7n 



Sr 



medium 



Sr 



(3.2) 



fine 



where pr* is some positive number. In practice, this ratio 
was taken to be 1.5 or 2. Each time the resolution is 
increased, the numerical solution must converge to the 
true solution. The numerical solution will have an error 
of the order of ( 5r*)^, then we can say that 

coarse = true + k{ Sr* f (3.3) 
'^medium = ^ true + K^^* / Pr'T (3-4) 
■^frne = true + k{Sr* / pl^f (3.5) 

where fc is a constant. From these relations it is easy to 
verify that 



fine 



medium 



medium 



(3.6) 



Figure 13 shows the differences (|3.6|) in absolute value. 
It corresponds to the instant t = 200M. By this time the 
initial pulse has bounced off the potential barrier and has 
reached a maximum amplitude of ~ 5 x 10^. A rough 
estimation indicates an error of 0.08% at the highest am- 
plitud. Fig. 21 shows the same phase than the previous 
one but &i t — lOOOAf. The error has increased now to 
^ 0.6%. A carefuU examination of Fig. 0] shows that the 
difference between the two lines is more notorious than 
in Fig. 13 That means that the convergence ratio is less 
than 4, but it is still consistent within a ratio of 3.95. 
Fig. [5l shows the same phase of the previous two graphs 
at different time steps. It is very clear that the relative 
error increases linearly with time. 

As for the convergence in the 9 direction, a minor prob- 
lem needs to be solved before computing the differences 
of the three numerical solutions. The problem is that the 
implementation of a staggered grid in 9 causes that the 
discrete set of values 9k were completely different when 
resolution is changed. This fact is illustrated in Fig. 
Therefore, in order to assess convergence, we used a sixth 
order Lagrange polynomial to interpolate the solution ob- 
tained with the medium and finer resolutions at the val- 
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FIG. 4: Convergence test in r* at 6 = n/2 and t = lOOOA/. 
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FIG. 5: Convergence test in r* sd, = ■k/2 at time intervals 
of lOOM. 



ues 6k of the coarse one. Having done that, Fig.[7|shows 
fourth order convergence for a fixed value of r* = 20M. 
The simulation parameters are still the same but this 
time Sr* = 0.25A/ and 69 = {tt/IG, 7r/24, 7r/36}, i.e. the 
ratio pg = 1.5. The amplitude decreases as the wave 
passes by. The graph shows different snapshot at inter- 
vals of 200M. Convergence is lost when round-off error 
is reached. 

Other test performed to check the validity of the nu- 
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FIG. 6: Distribution of the grid points 6k for different resolu- 
tions. 



FIG. 7: Convergence test in 9 at r* — 20M at time inter- 
vals of 200M. The continuous lines represent the fine minus 
the medium resolution (times 1.5'^) and the dashed one the 
medium minus the coarse one. 



merical solutions was evolving the initial data of a known 
analytic function. This function does not need to be a so- 
lution of the Teukolsky equation, it could be any smooth 
function in r* and 6, provided that the corresponding 
source term is added to the evolution equations. The idea 
is the following: let's call T the "Teukolsky operator" so 
that T($(t,r*,6')) = is the Teukolsky equation ifTTEIl . 
If we choose an arbitrary smooth function r* ,9), the 
result of applying the operator T will be 



Tmt,r*,9))^f{t,r*,9). 



(3.7) 



If we add the source term f{t,r*,9) to the evolution 
equations and give the initial data as $(0, r*,0) and 
9t$(O,r*,0) our code should reproduce the function 
In Fig. IHlwe show the result of such a test. We set the 
function $ as a Gaussian pulse (in r* and 9), traveling 
in the increasing direction of r*. Our code reproduce 
the analytic function with high accuracy. There is some 
damping in the amplitude of the pulse due to numeri- 
cal dispersion; however, this effect appears at very late 
times. 



C. Power-law tails 

The main application of this work is to accurately com- 
pute the power-law falloff in the gravitational perturba- 
tions evolution. In the following results the observation 
point is located at r* = 20M and 9 = Tr/2. The high- 
est resolution used was Sr* — 0.125M, 69 ~ tt/A8 and 
the lowest one was Sr* — IM, 66 = tt/S. Variation of 
the resolution in those intervals was done in order to ver- 
ify convergence, although the lowest resolution was not 
enough in cases where the 9 profile presents several os- 
cillations. Computationally, the determination of power- 
law tails is a challenging problem, because the amplitude 




of the wave decays exponentially during the quasinormal 
ringing and as an inverse power of time during the tail 
phase. This means that we are working with very small 
numbers that eventually reach round-off error, due to the 
finite precision of the computer processor. This posses 
some difficulties when we try to determine the exponent 
of the power-law. Recall that at very late times, for a 
finite value of r* (timelike infinity), the amplitude of the 
field goes as 

In principle, finding this exponent should not be a prob- 
lem since the field is proportional to a power of the time 
t. Thus, a simple power fitting of the form $ = At^^ 
(where A and fj, are constants) would be just enough. 
The problem with this procedure is that (|3.8|) is the 
very last stage in the evolution of the perturbation. In 
practice, we are not able to evolve the perturbations for 
such a long time, due to the finiteness of computational 
resources|33- So we analyze the field $ just from the 
moment the tail phase begins until the moment the solu- 
tion reaches round-off error. During this period, the field 
falloff is governed also by powers of time smaller than 
-(2/ -I- 3) Ml- 

$cxt-^-fO(i-(^+i)) (3.9) 

which means that the exponent of t reaches —{21 + 3) in 
an asymptotic way. Taking this fact in consideration, we 
compute the "local power index" defined as fiN = 
—tdt^/^, and use a linear fit (least-squares) such that 

fiN ^ fJ'+ — + (3.10) 

where the i?'s are constants. 

The least-squares fit yields good results when the lo- 
cal power index /xjy (that is computed taking numerical 
derivatives of <&) is taken in a large interval over which its 




FIG. 10: Power-law tails for m — 0, a = and different 
values of I. The evolution time is 1600M with a resolution 
Sr* = 0.25M and 30 = 0.098. 



oscillations are small. This behavior of at is illustrated 
in Fig. El For times longer than those shown in this plot, 
the oscillation of /zat becomes larger and the meaning of 
the local power index is lost. 

We notice that the larger the value of / the shorter the 
interval of validity for /ijv- This is due to the fact that 
for larger values of /, the power-law exponent is bigger 
(in absolute value), making the field to decay very fast 
reaching round-off error earlier. Fig. 1101 shows this be- 
havior for m = 0, a = and different values of Z. A 
closer examination of Figs. 1^ and 1101 reveals that the os- 
cillations in /i^v start some time before round-off error 
appears. The origin of such behavior is attributed to the 
accumulated numerical error that increases as the evolu- 
tion progresses. This oscillations are magnified in Fig. 
due to the numerical time derivatives of the field <f>. 

Table llVl shows the results a linear fit of the form (|3.10|) 
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I m predicted ^ ^ 0{t'^) ]rO(ir'~^r^) ^ 0{t''^) 



2 7 6.866 ± 0.004 7.00 ± 0.03 7.18 ± 0.01 

2 1 6.87 ± 0.01 7.01 ± 0.03 7.177 ± 0.003 

2 2 6.87 ± 0.01 7.04 ± 0.06 7.180 ± 0.005 

3 9 8.64 ± 0.04 8.9 ± 0.3 9.24 ± 0.02 
3 1 8.62 ± 0.04 8.9 ± 0.2 9.22 ± 0.02 
3 2 8.62 ± 0.02 8.9 ± 0.1 9.23 ± 0.01 

3 3 8.65 ± 0.04 9.1 ± 0.2 9.23 ± 0.02 

4 11 10.1 ± 0.7 11.3 ± 0.7 11.23 ± 0.04 
4 3 9.9 ± 0.1 10.7 ± 0.7 11.06 ± 0.04 
4 4 10.2 ± 0.1 11.5 ± 0.5 11.31 ± 0.03 



TABLE IV: Numerically computed power-law tails, using a 
least-squares fit to ^n- 
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TABLE V: Coefficients computed using the linear fitting. 



for the case a — Q and different values of I and m. The 
fourth column shows the value of the exponent /i using a 
fitting curve of the form /ijv = fi+Bi/t. The fifth column 
corresponds to a fitting curve fj,N = fJ, + Bi/t + B2/t'^ and 
the sixth column takes in consideration only a quadratic 
correction in t, i.e. /ijv = + B2 / ■ The uncertainty is 
the statistically computed error for the parameter /i in 
the fitting. This quantity is greater in the case of cor- 
rections 0{t~^ ,t~^) because there are three constants to 
be adjusted. The exponents are in agreement with the 
expected value 21 + 3. The agreement is better for ^ = 2 
than for I = 4 and for the 0{t~^ ,t~^) correction than 
for the 0(t^^) and 0{t~^) one. In two cases (/ — 4, 
m = 1,2), the duration of the tail was not enough to 
determine the exponent. The interval chosen to make 
the curve fitting was the largest one that starts after the 
quasinormal ringing and ends before the amplitude of the 
UN oscillations got too high in a way that it could bias the 
result. In tabled we show the values of the coefficients 
Bi and B2 for the three different fitting functions. It is 
very interesting to notice that the weight of the term 
is greater than that of by a factor of at least 200. This 
result supports, to some extend, the model proposed by 
Poisson [l^l for the radiative falloff of a scalar field in a 
stationary, asymptotically flat and weakly curved space- 
time. He shows that the first correction to the power-law 
tail is of order t^^. Our observation point is located at 
r* = 20Af and it is not far enough from the black hole 
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TABLE VL Power law tails, scalar case I = 0. Initial data: 
outgoing pulse. (Observer position in M units) 
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FIG. 11: Power law tails, scalar case Z = 0,m = 0. Observer's 
position is in M units. 

to say that it is in the asymptotically flat region. 

To find out the behavior in the asymptotically flat re- 
gion we did similar evolutions for the scalar case s = 
0,; = and TO = 0. Tables IVll and IVlIl show the power 
law tail computed at different observation points in the 
equatorial plane. In the first one the initial data was an 
outgoing Gaussian pulse centered at r* — lOOM. In the 
second one, initial data has the same initial shape but 
the pulse has zero velocity, i.e. it is time symmetric. We 
see that the power law tail is 4 instead of 3 for time sym- 
metric initial data. The fact that the fitting to a function 
of order yields tails closer to the predicted value for 
distant observer supports Poisson's formula. Figure 1111 
show the exponential fall off at different distances from 
the black hole as a function of time. All them approach 
asymptotically to the theoretically known value. 

Figure IT^ shows the evolution for the I = 2 multipole 
for different values of to. In these cases there is no pres- 
ence of round-off error because it has the slowest decay 
rate, $ oc In Fig. ^| we see the same situation as 
above but with 1 — 3. The power-law has a behavior 
$ (X t~^. Round-off error appears at $ ~ 10~^^. Fi- 
nally in Fig. ll4l round-off error appears approximately at 
the same value of $ as in the previous case. We notice 
that the quasinormal ringing is practically inexistent. In 
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o{t-\t-^) 


20 


3.959 


4.005 


4.001 


50 


3.950 


4.005 


4.002 


150 


3.859 


3.999 


4.019 


250 


3.652 


3.981 


4.079 


350 


3.601 


3.979 


4.084 


450 


3.296 


3.941 


4.219 


550 


3.419 


3.961 


4.143 


650 


3.128 


3.921 


4.291 



TABLE VII: Power law tails, scalar case I = 0. Initial data: 
zero velocity pulse. (Observer position in M units) 




FIG. 12: Power-law tail Z = 2, m = 0, 1, 2, a = 0. 



all the runs we use outgoing initial data as prescribed in 

So far we have been considering power-law tails for the 
case of a non-rotating black hole (a = 0). Fig. 1151 shows 
the power-law tails for the case in which a — 0.5. We 
can see that for these values of I the duration of the tail 
is too small. Doing a nonlinear fitting for the case I = 2, 
we found a tail of —7.0011. For the other cases, this very 
short tail phase is not enough to tell with certainty the 
power-law exponent; besides, the tail has still some small 
oscillation in that time interval. To verify the effects of 
a Kerr spacetime in the evolution of the gravitational 
perturbations, the frequencies of the quasinormal ringing 
are useful; as shown in the next section. 



D. Quasinormal modes 

We compute the quasinormal modes frequency for the 
cases shown in Fig. ^1 which correspond to the Kerr 
spacetime ftable IVIIl| . The case I = 2 agrees with the 
known frequencies |23l with an error less than 1%. 
The cases I = 3 and I = 4 have frequencies similar to 
the I = 2 multipole. We also compute these frequencies 
for the case of a Schwarzschild spacetime ftable llX|) . In 





FIG. 14: Power-law tail « = 4, m = 0, 1, 3, 4, a = 0. 




FIG. 15: Power-law tails for a = 0.5. 
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FIG. 16: Power-law tail I = 3, m = 0, a = 0, 0.5, 0.9. 
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FIG. 17: Power-law tail Z = 4, m = 0, a = 0, 0.5, 0.9. 



this case, the frequency values agree with the predicted 
ones IMEEEl] within a 0.1 to 1.4% error. In the cases 
for I = 4, the numerical evolution was not able to render 
quasinormal ringing. For a ^ 0, angular mode conversion 
to the lowest allowed multipole is present. The case I — 
3, TO = presents an irregular oscillation and the real 
part of the frequency and cannot be trusted. That is the 
reason why it is not shown here. 

Finally, Figs . 1 1 61 and 1 1 71 show the evolution for m = 



/ 


m 


computed a M 


predicted aM 


2 





0.384 + 0.0875i 


0.3833 -1- 0.08707i 


2 


1 


0.341 + 0.0805i 


0.4206 + 0.08617i 


3 





NA -HO.OSOOi 


0.61212 + 0.09077i 


3 


1 


0.339 + 0.0803i 


0.65060 + 0.0900i 


4 





0.382 + 0.0860i 


NA 


4 


1 


0.341 + 0.0797i 


NA 



TABLE VIII: Quasinormal mode frequencies for a = 0.5 (NA 
= not available). 



1 m a M predicted ctM 

2 0.373 -I- 0.0875i 0.3736715 + 0.0889625i 
2 1 0.375 -I- 0.0869i 0.3736715 + 0.0889625i 

2 2 0.376 + 0.0877i 0.3736715 + 0.0889625i 

3 0.60H- 0.0903i 0.5994435 0.092703i 
3 1 0.600 + 0.0902i 0.5994435 0.092703i 
3 2 0.605 -I- 0.090H 0.5994435 -^- 0.092703i 
3 3 0.598 + 0.0933i 0.5994435 -f 0.092703i 



TABLE IX: Quasinormal mode frequencies for o = 0. 

of / = 3 and / = 4 respectively, for different values of a. 
We can see roughly that the power-law tail is the same 
for each case and that multipole conversion is present. 

E. Fourth order versus second order 

We could say that the advantage of using a fourth order 
convergent code is that we can achieve the same results of 
the second order one with less resolution. In other words, 
the same degree of accuracy can be obtained with both 
approaches in the same computational domain but the 
second order one will need more points. Quantitatively, 
we can compare the error in the solution for both cases. 
This error is the difference between the true solution and 
the numerical solution for a given resolution. If we denote 
this quantity by en, where h is the grid spacing then we 
have: 

eh = fc„/i" (3.11) 

where /c is a constant and n is the order of accuracy. If we 
equate the errors for n — 2 and n = 4, the relationship 
between resolutions is 

h2 = '^hl (3.12) 

If in the fourth order method, ^,4 = 0.1 then the equiva- 
lent resolution in the second order one is approximately 
10 times bigger! i.e. h2 — 0.01. This implies that in a 
one dimensional problem the number of points is also 10 
times bigger. If a two dimensional problem is considered 
then the second order grid should contain 100 times more 
points than that of the fourth order method, to get the 
same error in the solution. 

A feature of the finite differences methods is that ac- 
cording to (|3.11|) , if we increase resolution by some factor 
c the error is reduced by a factor c". Thus each time we 
double resolution, the error decreases by a factor of 4 in 
the second order method and by 16 in the fourth order 
method. 

The above considerations put the fourth order method 
in a better position than the second order one but, as we 
said in the introduction, the price we have to pay is run- 
ning time. Given a resolution h, the fourth order method 
will find a more accurate solution than the second order 
one. The time that the fourth order method will take 
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RAM memory [Mb] running time [hrs] 
Nr X Ne 2nd order 4th order 2nd order* 4th order 



1000 X 8 


19 


22 


0.06 


0.25 


2000 X 16 


64 


71 


0.48 


2 


4000 X 32 


236 


249 


3.73 


16 



TABLE X: RAM memory used grid different grid sizes for 
both second and fourth order methods. *These times are for 
equal grid sizes. For equivalent resolutions, the running time 
for the 2nd order method is approximately 14 times larger 
than that of the 4th order method. 

will be longer because there are much more calculations 
to be done. In order to determine if the gain in a smaller 
grid is greater than the loss in running time, we did some 
numerical experiments. The running time tmn is given 
approximately by 

trun-^(^^^ to, (3.13) 

where io is the running time at resolutions Sr^ and S9q. 
Fixing 6rQ — 1 and dd^ — vr/S, ~ 16.3 min in the case 
of the fourth order method; and to — 3.5 min for the 
second order one (we used a Pentium 4 CPU 2.4GHz). 
These times correspond to the gravitational case 1 = 2, 
m = in an interval — 50 < r* < 950M being the simu- 
lation time 1600 M. In this simulation, the fourth order 
method gave very good results whereas the second or- 
der one becomes unstable around t = 300M. In order to 
obtain the same result (power-law tail) than the fourth 
order one, it was necessary the increase the resolution 
four times in both directions. This implies that now, the 
nmning time for the second order method is tmn ^ 230 
min. This time is 14 times larger than that correspond- 
ing to the fourth order method. So we definitely have a 
gain in speed, when the errors in the numerical solutions 
(for both methods) are kept equal. 

As to the RAM memory, table El gives information 
about the amount of memory used in function of the grid 
size {Nr X Ng). These values depend on the coding de- 
tails of the algorithm. In this kind of problem, computer 
memory is not a crucial factor as it is the speed. 



IV. CONCLUSIONS 

The main goal of this research has been to imple- 
ment a stable fourth order accurate method to numer- 
ically integrate the Teukolsky equation in the time do- 
main. In order to verify and evaluate the efficacy of 
our fourth order method, we have reproduced the main 
known results, i.e. power-law tails and quasinormal ring- 
ing. Power-law tails is a subject in which there are 
still unsolved questions concerning the late time behavior 
of the perturbations and its dependence on the coordi- 
nates and the initial data 0. We addressed some of 
these issues here and others are left for future research. 



Among the formers, we studied the case of an initial pulse 
with an initial configuration such that ^^l*^" ^ and 
3t<&£|*"" = 0. We have been able in this case that the late 
time behavior agrees well with the predicted Tl'l decay 
^ l/i^^+'*. We have also confirmed numerically that for 
an observer located far away from the hole {robs > 20Af) 
we see the predicted lOj correction to the power law as 

The fourth order method implemented in this work 
has shown to be convergent and stable in all runs we did. 
The ansatz proposed by Krivan et. al. [3 effectively 
removes the growing in time of the field, a feature that is 
expected from the asymptotic behavior of the solutions. 
When carrying out our calculation, we notice that the 
coeSicients of the Teukolsky equation reported in 0| do 
not correspond to those we got when the ansatz proposed 
is used. Instead, they correspond to an ansatz without 
the factor (see eq. H2.3|l '). We also verify that the for- 
mulation of Teukolsky equation as a system of first order 
differential equations is a powerful technique because it 
allows to express the first time derivative in function of 
the spatial derivatives. The resulting system of equa- 
tions has the form of the advection equation. That was 
the key idea that allowed us to implement a fourth or- 
der method to solve the equation: expand the solution 
in power series of time keeping terms up to fourth order 
and then use the differential equation to substitute the 
time derivatives by spatial derivatives. Using this pro- 
cedure, a higher order method could be developed in a 
straightforward manner. The stability analysis provided 
a limiting value of the Courant of 1.5 below which we can 
perform stable evolutions. 

The time evolutions carried out with the fourth order 
method yielded accurate results even when relatively low 
resolutions were used. The lowest resolution was Sr* = 1, 
66 = tt/8 for the I = 2, m = case. The highest reso- 
lution used was Sr* = 0.125 and 66 = 7r/48, for the 
I = 4 multipole. In finite difference methods the res- 
olution is chosen in such a way that the details in the 
profiles of the functions involved in the calculations can 
be accurately approximated. For higher values of Z, the 9- 
dependence has more oscillations in its domain therefore 
more points are needed to find a reliable solution. The 
more oscillations or narrow peaks a function has the more 
resolution we need. This comes from the fact that those 
functions have derivatives whose values oscillate rapidly 
and higher derivatives vary even faster. If the resolution 
is not good enough the effect is a "numerical mode mix- 
ing" , that has nothing to do with the physical model, but 
with the numerical aspects of the implementation. This 
mode mixing acts like if we were evolving the waves in 
a Kerr background, where physics tells us that angular 
mode mixing is expected. This unwanted effect cannot be 
easily detected when a 7^ in the simulations. Therefore 
it is absolutely necessary to verify that this effect is not 
present when we evolve in the Schwarzschild background, 
where the physical mode mixing does not happen. This 
may explain some discrepancy on the computed power- 
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law exponent decay appeared in the literature. 

For a ^ table IVTTT1 shows clearly that for ^ = 3,4, .. 
the mode mixing acts bringing down the quasinormal fre- 
quencies close to that of the 1 = 2. The same effect is 
observed in the tails power, although it is more difficult 
to prove (See Fig.[TB|l. 

Boundary conditions were not an issue of concern in 
this research. The radial inner boundary was not a prob- 
lem because the field decays exponentially near that re- 
gion. The radial outer boundary always reflects part of 
the wave. The immediate solution is to push this bound- 
ary far away such that this reflection does not interfere 
in the region of interest. A refinement of the boundary 
conditions will be needed when a large computational 
domain can not be used in favor of higher resolutions. 

We have a reliable computational tool to explore sev- 
eral interesting problems concerning first order pertur- 
bations. We expect to shed some light on problems like 
the late time behavior of gravitational and scalar fields 



in a Kerr background. This is a problem that has been 
studied both analytically and numerically in recent years 
and there are still questions to be answered. Further re- 
search could include the problem of the orbiting particle 
around a black hole and the close limit a ppr oximat ion in 
the problem of two colliding black holes |28l |. 



Acknowledgments 

We are very grateful to E. Poisson for useful discussions 
on power-law tails, and to K. Kokkotas for making avail- 
able unpublished values of the Quasinormal frequencies 
for rotating black holes. We also gratefully acknowledge 
the support of the NASA Center for Gravitational Wave 
Astronomy at The University of Texas at Brownsville 
(NAG5-13396) and also to NSF for financial support from 
grants PHY-0140326 and PHY-0354867. 



I 

APPENDIX A: FOURTH ORDER ACCURATE DERIVATIVES 



This are the formulas to compute fourth order accurate derivatives using finite differences. Those which are centered 
have always less truncation error than the corresponding off-centered ones. In all of them, h is the size of the step 
and a; < ^ < a; + ft,. We neglect term of powers higher than four. 



• First derivative 
— centered: 



off-centered (1 point): 



^ -3u,-_i-10^.,- + 18.,+i-6^.,-+2 + ^,+3 1^4^(5) (A2) 
^ ' 12 h 20 ^ ' 



off-centered (2 points): 



u'ix) = -25 u, + 48 u,+, 36 + 16 ^.,+3 - 3 u,+, 1^4^(5) (^3) 

12 h 5 



• Second derivative 
— centered: 



u'\x) = -".--2 + 16^,-i-30^,- + 16.,-+i-z.,-+, _ 1^,4^(6) (A4) 

1 2 /i 90 



— off-centered (1 point): 
u"{x) 



„ 10uj_i — 15uj — 4wj+i + 14uj+2 — 6uj_|_3 + 



12 

— h^u^'HO (A5) 
180 ^ ' 



— off-centered (2 points): 

45 Uj — 154 Uj+i + 214 Mj+2 — 156 Wj+a + 61 itj+4 — 10 



u {x) 



12ft2 

_H!/jV6)(C) (A6) 
180 ^ ' 



Third derivative 
— centered: 



_ ~ 8 Uj_2 + 13 - 13 Uj+1 + 8 - Wj_|_3 



— off-centered (1 point): 

v!'\x) 



1,1 , . _ —u.j-2 — 8 Uj-i + 35 u-j — 48 Uj+i + 29 — 8 uj+s + u^+i 



8h^ 



off-centered (2 points): 

— 15 Uj-i + 56 Uj — 83 Uj+i + 64 — 29 Wj+s + 8 — Wj+s 



u'"(a;) 



8/i3 



off-centered (3 points): 

-49 Uj + 232 Uj+i - 461 Uj+2 + 496 uj+s - 307 104 Uj+5 



u"'{x) 



8/i3 



Fourth derivative 
— centered: 



(4). ^ _ —Uj-3 + 12 Uj-2 — 39 Uj-1 + 56 Uj — 39 Mj+i -|- 12 Uj+2 - Uj+3 
^"^^ ~ 67? 
7 



240 



off-centered (1 point): 



(4), . _ 4 Uj-2 - 11 Mj-i + 31 Uj+i - 44 Mj+2 -h 27 Uj+3 — 8 Mj+4 + Uj+5 
^""^ ~ 67? 



off-centered (2 point): 



21 Uj-i — 112 Uj + 255 Uj+i — 324 Uj+2 + 251 Uj+3 — 120 Uj+4 
= 67? 

33Mj+5-4Uj-+6 _ 127 4 (8), . 

+ 6/i4 240 



off-centered (3 points): 



(•4-) 56 Uj — 333 Mj+i + 852 w.j+2 — 1219 Mj+3 -|- 1056 Mj+4 — 555 w^+s 
= 67? 



164mj+6-21mj+7 _ 967^4^ (8),^, 
^ 6/i4 240 
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